/* Colluding Against Workers: Main model: bootstrapping
	Delabastita & Rubens 
------------------------------------------*/

clear
set more off
cd "$station"

 
set matsize 1000 
local B = $bs		 // number of bootstrap iterations
 
di `B'
  
set matsize 1000	
set more off
clear  

** Matrices to store estimates
 
forvalues m = 1(1)1{
foreach mod in   "bb" "ols" "rts" {
foreach coef in "thl_`mod'" "thm_`mod'"  "scale_`mod'" "bl_`mod'" "rho_`mod'" /* "bd_`mod'"*/  "bk_`mod'" "blk_`mod'" "bmk_`mod'"   "blm_`mod'"  "bm_`mod'" {
matrix _`coef'`m' = J(`B', 1, 0)
matrix m_`coef'`m' = J(`B', 1, 0)
}
foreach mod in   "bb" "ols" "rts" {
foreach coef in "md" "mu"  "conv"  "mup" "mdp"  {
matrix _`coef'_`mod'`m' = J(`B', 1, 0)
matrix m_`coef'_`mod'`m' = J(`B', 1, 0)
matrix a_`coef'_`mod'`m' = J(`B', 1, 0)
}
}
}

foreach n in "a" "b" "c" {
forvalues m = 1/3 {
matrix _mdcorr`m'`n' = J(`B', 1, 0)
}
}

forvalues m = 1/4 {
matrix _kappa_dq`m' = J(`B', 1, 0)
matrix m_kappa_dq`m' = J(`B', 1, 0)
matrix _kappa_rts_dq`m' = J(`B', 1, 0)
matrix m_kappa_rts_dq`m' = J(`B', 1, 0)
}
  
forvalues m = 1/2 {
foreach coef in   "th_rail`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'"  "th_car`m'"   ///
 "th_lib`m'" "th_soc`m'" "th_com`m'" "th_oth`m'" "th_libp`m'" "th_socp`m'" "th_comp`m'" "th_othp`m'"   { 
matrix _`coef' = J(`B', 1, 0)
matrix m_`coef' = J(`B', 1, 0)
} 
}

forvalues m = 1/4 {
foreach coef in     "th_union`m'" ///
   { 
matrix _`coef' = J(`B', 1, 0)
matrix m_`coef' = J(`B', 1, 0)
} 
}

  
foreach coef in   ///
 "th_yr1" "th_yr2" "th_yr3" "th_yr4" "th_yr5" "th_yr6"  "th_lyr"    { 
matrix _`coef' = J(`B', 1, 0)
matrix m_`coef' = J(`B', 1, 0)
} 
}
  
use ./data/temp/data_pf, clear 
di _N
 
** Start bootstrapping iteration

forvalues b = 1 /`B' {
preserve
set seed `b'
bsample   ,  cluster(mineid) idcluster(idb)
drop mineid
rename idb mineid   

xtset mineid yr
 
/* A. Production function */ 

** Cobb-Douglas   

xtset mineid yr
 
* OLS

gen lsegers = log(segers_nomwage_agr / cpi)
reg lq lemp lwm  lk  

mat coef = e(b)
mat list coef
gen rho_ols1 = .
gen bl_ols1 = coef[1,1]
gen bm_ols1 = coef[1,2]
gen bk_ols1 = coef[1,4]
gen bc_ols1 = coef[1,5]

gen thl_ols1 = bl_ols1  
gen thm_ols1 = bm_ols1  
gen scale_ols1 = bl_ols1 + bm_ols1 + bk_ols1 

* GMM

foreach var of varlist segers_nomwage_* {
	gen l`var' = log(`var')
}

sum winv*
gen lwinv = log(winv+1)
gen lud = log(u/(s+u))

gen lemp2 = lemp^2
gen lwm2 = lwm^2
gen lk2 = lk^2

gen lemp3 = lemp^3
gen lwm3 = lwm^3
gen lk3 = lk^3

gen lsegers2 = lsegers^2

gen lwinv2 = lwinv^2

 gen lsegers_ag = log(segers_nomwage_agr/cpi)
 gen lsegers_ser = log(segers_nomwage_ser/cpi)
 gen lsegers_ind = log(segers_nomwage_ind/cpi)
  
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )      , ///
inst(  L.lemp  L.lwm  L.lsegers L.lk lk  )  conv_maxiter(50)  

  
display as text "iteration = "   `b'

gen conv = e(converged)

mat coef = e(b)
mat list coef
gen rho_bb1 = coef[1,1]
gen bl_bb1 = coef[1,2]
gen bm_bb1 = coef[1,3]
gen bk_bb1 = coef[1,4]
gen bc_bb1 = coef[1,5]
 gen conv1 = e(converged)  
gen thl_bb1 = bl_bb1  
gen thm_bb1 = bm_bb1  

foreach var of varlist *bb* {
	replace `var' = . if conv==0
}

gen scale_bb1 = thl_bb1 + thm_bb1 + bk_bb1
   
  ** Fix returns to scale to nu = 1.05

gen nu = 1.05
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -(nu-{bl}-{bm})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk L.lsegers  )  conv_maxiter(50) 
 
local N_rts = e(N)
 
mat coef = e(b)
mat list coef
gen rho_rts1 = coef[1,1]
gen bl_rts1 = coef[1,2]
gen bm_rts1 = coef[1,3]
gen bk_rts1 = 1-bl_rts -bm_rts 
gen bc_rts1 = coef[1,4]
gen blm_rts1  = .
gen blk_rts1 = .

gen thl_rts1  = bl_rts  
gen thm_rts1  = bm_rts 
gen scale_rts1 = nu 

 estat overid 
 local j  = string(r(J)  ) 		 
local j_rts1 = substr("`j'" ,1,4) 

local jp  = string(r(J_p)  ) 		 
local jp_rts1 = substr("`jp'" ,1,4)  
      
/* B. Markdowns and markups
---------------------------*/

gen am = wm/rev
gen al = wl/rev
 
foreach mod in "bb" "ols" "rts" {
gen mu_`mod'1 = bm_`mod'1/am
gen md_`mod'1 = ((bl_`mod'1)/al)/mu_`mod'1
gen mp_`mod'1 = 1 + (mu_`mod'1-1)+(md_`mod'1-1)
gen mdp_`mod'1 = (md_`mod'1 - 1)/md_`mod'1
gen mup_`mod'1 = (mu_`mod'1 - 1)/mu_`mod'1
}

 
bys yr: egen agemp = sum(emp)
gen sempag = emp/agemp

 
** Average markdowns per year
 
foreach var in "mu" "md"  "mdp" "mup" {
foreach mod in "bb" "rts" "ols" {
forvalues m = 1/1 {
bys yr: egen  m`var'_`mod'`m' = median(`var'_`mod'`m')
bys yr: egen  av`var'_`mod'`m' = mean(`var'_`mod'`m')
bys yr: gen s`var'_`mod'`m' = `var'_`mod'`m'*sempag
bys yr: egen  ag`var'_`mod'`m' = sum(s`var'_`mod'`m')
}
}
}

sum mmd_rts1 , d 

  
  
forvalues m = 1(1)1 {
foreach mod in   "bb" "ols" "rts" {
foreach var in  "thm_`mod'" "thl_`mod'" "bk_`mod'" "bl_`mod'" "rho_`mod'"   "bm_`mod'"  "scale_`mod'" {
 egen av`var'`m' = mean(`var'`m')
}
}
}
 
 


** Concentration
 
xtset mineid yr
drop open
gen open= q>0
replace open = . if q==.
sum yr
   
 forvalues m = 1(1) 4 {
bys dq`m'id yr: egen n_dq`m' = sum(open)
tab n_dq`m' , gen(nm_dq`m'_)
}

gen urb = dq2=="40" | dq2=="50"
bys mineid: egen urban = max(urb) 
  
  foreach var of varlist md* mu* {
gen l`var' = log(`var')
}


capture: reg lmd_bb1 nm_dq3_1 nm_dq3_2 nm_dq3_3 drail dtram urban  i.yr, vce(cluster dq3id )

 
 capture: gen th_rail1 = _b[drail]
capture: gen th_tram1 = _b[dtram]
capture: gen th_nm11 = _b[nm_dq3_1]
capture: gen th_nm21 = _b[nm_dq3_2]
capture: gen th_nm31 = _b[nm_dq3_3]
capture: gen th_urban1 = _b[urban]

capture: xtreg lmd_bb1 nm_dq3_1 nm_dq3_2 nm_dq3_3 drail dtram urban   i.yr, fe r
capture: gen th_rail2 = _b[drail]
capture: gen th_tram2 = _b[dtram]
capture: gen th_nm12 = _b[nm_dq3_1]
capture: gen th_nm22 = _b[nm_dq3_2]
capture: gen th_nm32 = _b[nm_dq3_3]
capture: gen th_urban2 = _b[urban]



* Variance decomposition

bys mineid: egen e_car=max(dcar)
 /*
gen luw = log(wl/emp)
reg luw i.yr  
gen r2_11 = e(r2)
reg luw i.yr i.dq4id  
gen r2_21 = e(r2)
reg luw i.yr i.dq4id  
gen r2_31 = e(r2)
areg luw i.yr i.dq4id  , absorb(mineid)
gen r2_41 = e(r2)

capture: reg lmd_bb1 i.yr
capture: gen r2_12 = e(r2)
capture: reg lmd_bb1 i.yr i.dq4id  
capture: gen r2_22 = e(r2)
capture: reg lmd_bb1 i.yr i.dq4id  
capture: gen r2_32 = e(r2)
capture: areg lmd_bb1 i.yr i.dq4id , absorb(mineid)
capture: gen r2_42 = e(r2)
*/

/* Labor supply model */
  
* Cournot
   bys mineid: egen e_union=max(dunion)
save ./data/temp/data_cournot, replace

forvalues m = 1/4 {
use  ./data/temp/data_cournot, clear
collapse(sum) emp wl drail dtram e_union dcar e_car cpi   degreve*, by(dq`m'id   yr)
gen p98 = yr>=1898
gen dcar_p98 = dcar*p98
gen e_car_p98 = e_car*p98
gen demshock = yr>= 1871 & yr<= 1875
gen lemp = log(emp)
gen luw = log(wl/emp)
gen lqimp = log(degreve_imp_q)
gen pimp = (degreve_imp_p/cpi)
gen lpimp = log(pimp)

reg lqimp lpimp
predict demshock2, resid

gen ptrans =  ((degreve_trans_pq/degreve_trans_q)/cpi)
gen lptrans = log(ptrans)
ivregress 2sls luw (lemp =  demshock  e_car_p98  ) e_car p98  
global bw_dq`m' = _b[lemp]
}

use  ./data/temp/data_cournot, clear

gen bw_dq1 =  $bw_dq1
gen bw_dq2 =  $bw_dq2
gen bw_dq3 =  $bw_dq3
gen bw_dq4 =  $bw_dq4
   
 
forvalues m = 1/4 {
bys dq`m'id yr: egen emp_dq`m' = sum(emp)
gen semp_dq`m' = emp/emp_dq`m'
gen psic_dq`m' = 1+bw_dq`m'*semp_dq`m'
gen psim_dq`m' = 1+bw_dq`m'
gen kappa_dq`m' = (md_bb1- psic_dq`m'    )/( psim_dq`m' - psic_dq`m') 
gen kappa_rts_dq`m' = (md_rts1- psic_dq`m'    )/( psim_dq`m' - psic_dq`m') 
foreach var in "psic" "psim" "kappa" "kappa_rts"{
bys yr: egen av`var'_dq`m' = mean(`var'_dq`m')
bys yr: egen med`var'_dq`m' = median(`var'_dq`m')
}
}



** Collusion: cross-sectional

gen lkappa_dq4 = log(kappa_dq4+1)
capture: reg lmd_bb1 e_union dcar yr , r

capture: gen th_union1 = _b[e_union]
capture: gen th_car1 = _b[dcar]
 
capture: reg lmd_bb1 _dq4 e_union dcar i.yr , r
capture: gen th_union2 = _b[e_union]
capture: gen th_car2 = _b[dcar]


** Collusion and employer associations: comparison pre- and post-cartel period:

capture:  reg lmd_bb1  e_union  i.yr i.dq1id if yr<1898 & dcar~=., r 
capture: gen th_union3 = _b[e_union]

capture: reg lmd_bb1  e_union  i.yr i.dq1id if yr>=1898 & dcar~=., r
capture: gen th_union4 = _b[e_union]


** Markdown-size correlation

gen lsemp_dq4 = log(semp_dq4)
egen mid = group(dq4id yr)

* (a) all firms

capture: reg lmd_bb1 lsemp_dq4  yr
capture: gen mdcorr1a = _b[lsemp_dq4]
capture:  areg lmd_bb1 lsemp_dq4 yr, absorb(dq4id)
capture: gen mdcorr1b = _b[lsemp_dq4]
capture: areg lmd_bb1 lsemp_dq4 yr, absorb(mid)
capture: gen mdcorr1c = _b[lsemp_dq4]

* (b) non-cartel firms

capture: reg lmd_bb1 lsemp_dq4  yr if dcar==0
capture: gen mdcorr2a = _b[lsemp_dq4]
capture: areg lmd_bb1 lsemp_dq4 yr if dcar==0, absorb(dq4id)
capture: gen mdcorr2b = _b[lsemp_dq4]
capture: areg lmd_bb1 lsemp_dq4 yr if dcar==0, absorb(mid)
capture: gen mdcorr2c = _b[lsemp_dq4]

* (c) cartel firms

capture: reg lmd_bb1 lsemp_dq4  yr if dcar==1
capture: gen mdcorr3a = _b[lsemp_dq4]
capture: areg lmd_bb1 lsemp_dq4 yr if dcar==1, absorb(dq4id)
capture: gen mdcorr3b = _b[lsemp_dq4]
capture: areg lmd_bb1 lsemp_dq4 yr if dcar==1, absorb(mid)
capture: gen mdcorr3c = _b[lsemp_dq4]


 ** Collusion: changes over time

** Time trend 

capture: reg lmd_bb1  yr ,  r
capture: gen th_lyr = _b[yr] 

gen dyr1 = yr<1865 & yr>=1855
gen dyr2 = yr<1875 & yr>=1865
gen dyr3 = yr<1885 & yr>=1875
gen dyr4 = yr<1895 & yr>=1885
gen dyr5 = yr<1905 & yr>=1895
gen dyr6 = yr<1915 & yr>=1905

capture: reg lmd_bb1  dyr1 dyr2 dyr3 dyr4 dyr5 dyr6
capture: gen th_yr1 = _b[dyr1] 
capture: gen th_yr2 = _b[dyr2] 
capture: gen th_yr3 = _b[dyr3] 
capture: gen th_yr4 = _b[dyr4] 
capture: gen th_yr5 = _b[dyr5] 
capture: gen th_yr6 = _b[dyr6] 
 
 
 
sum medpsic* medkappa*
 
 
** Store estimates
 
forvalues m = 1/2 {
 foreach coef in       "th_rail`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'" "th_car`m'"   {  
 capture:   egen av`coef' = mean(`coef')
}  
}


forvalues m = 1/4 {
 foreach coef in       "th_union`m'"{  
 capture:   egen av`coef' = mean(`coef')
}  
}

foreach coef in   "th_yr1" "th_yr2" "th_yr3" "th_yr4" "th_yr5" "th_yr6"  "th_lyr"     {
 capture:   egen av`coef' = mean(`coef')
}  

keep in 1/`B'
 
* Time trend

forvalues m = 1(1)1 {
	foreach mod in "bb" "rts" "ols" {
		foreach var in    "md"  "mu" "mup" "mdp"  {
 mkmat av`var'_`mod'`m', matrix(av`var'_`mod'`m')
  mkmat m`var'_`mod'`m', matrix(m`var'_`mod'`m')
  mkmat ag`var'_`mod'`m', matrix(ag`var'_`mod'`m')
 matrix _`var'_`mod'`m'[`b',1] = av`var'_`mod'`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
  matrix m_`var'_`mod'`m'[`b',1] = m`var'_`mod'`m'[1,1]		 
  matrix a_`var'_`mod'`m'[`b',1] = ag`var'_`mod'`m'[1,1]		 
}  
}
}




foreach n in "a" "b" "c" {
forvalues m = 1/3 {
 capture: matrix _mdcorr`m'`n' = J(`B', 1, 0)
 capture: mkmat mdcorr`m'`n', matrix(mdcorr`m'`n')
 capture: matrix _mdcorr`m'`n'[`b',1] = mdcorr`m'`n'[1,1]			// store average markdown for this bootstrap iteration in matrix
}
}

forvalues m = 1/4 {
 mkmat avkappa_dq`m', matrix(avkappa_dq`m')
 mkmat medkappa_dq`m', matrix(medkappa_dq`m')
 mkmat avkappa_rts_dq`m', matrix(avkappa_rts_dq`m')
 mkmat medkappa_rts_dq`m', matrix(medkappa_rts_dq`m')
 }

   
forvalues m = 1(1)1 {
foreach mod in   "bb" "ols" "rts" {
foreach var in "thl_`mod'"  "thm_`mod'" "bk_`mod'" "bl_`mod'"    "rho_`mod'"  "bm_`mod'" "scale_`mod'"   {
  mkmat av`var'`m', matrix(av`var'`m')
 matrix _`var'`m'[`b',1] = av`var'`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
}
}
}

forvalues m = 1(1)4 {
foreach var in "kappa"  "kappa_rts" {
  mkmat av`var'_dq`m', matrix(av`var'_dq`m')
  mkmat med`var'_dq`m', matrix(med`var'_dq`m')
 matrix _`var'_dq`m'[`b',1] = av`var'_dq`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
 matrix m_`var'_dq`m'[`b',1] = med`var'_dq`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
}
}


forvalues m = 1/2 {
foreach coef in   "th_rail`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'" "th_car`m'"   ///
 { 
 capture:    mkmat av`coef', matrix(av`coef')
 capture:  matrix _`coef'[`b',1] = av`coef'[1,1]			// store average markdown for this bootstrap iteration in matrix
} 
}


forvalues m = 1/4 {
foreach coef in     "th_union`m'"  { 
 capture:    mkmat av`coef', matrix(av`coef')
 capture:  matrix _`coef'[`b',1] = av`coef'[1,1]			// store average markdown for this bootstrap iteration in matrix
} 
}
 
foreach coef in    "th_yr1" "th_yr2" "th_yr3" "th_yr4" "th_yr5" "th_yr6"  "th_lyr"       {
 capture:    mkmat av`coef', matrix(av`coef')
 capture:  matrix _`coef'[`b',1] = av`coef'[1,1]			// store average markdown for this bootstrap iteration in matrix
}  

 
 
restore			// end bootstrapping loop here
}
 
 
 
foreach n in "a" "b" "c" {
forvalues m = 1/3 {
svmat _mdcorr`m'`n' 
}
}
 
  
 
**********

forvalues m = 1(1)1 {
foreach mod in "bb" "ols" "rts" {	
foreach var in "md"   "mu"  "mdp" "mup"  {
svmat _`var'_`mod'`m'   
svmat m_`var'_`mod'`m' 
svmat a_`var'_`mod'`m' 
}
}
}



forvalues m = 1(1)4 {
foreach var in "kappa"  "kappa_rts"  {
svmat _`var'_dq`m'   
svmat m_`var'_dq`m'   
 }
}
 

forvalues m = 1(1)1 {
foreach mod in   "bb" "ols" "rts" {
foreach var in "thl_`mod'"  "thm_`mod'" "bk_`mod'" "bl_`mod'" "bm_`mod'" "blk_`mod'" "bmk_`mod'"   "rho_`mod'" /*"bd_`mod'"*/   "scale_`mod'" /*"blloc_`mod'" "blcut_`mod'"*/ {
svmat _`var'`m'   
}
}
}
 
forvalues m = 1/2 {
foreach coef in   "th_rail`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'" "th_car`m'"  ///
 "th_lib`m'" "th_soc`m'" "th_com`m'" "th_oth`m'" "th_libp`m'" "th_socp`m'" "th_comp`m'" "th_othp`m'"    { 
/*capture:*/ svmat _`coef' 		// store average markdown for this bootstrap iteration in matrix
} 
}
 forvalues m = 1/4 {
foreach coef in    "th_union`m'"    { 
/*capture:*/ svmat _`coef' 		// store average markdown for this bootstrap iteration in matrix
} 
 
}
 
foreach coef in   "th_yr1" "th_yr2" "th_yr3" "th_yr4" "th_yr5" "th_yr6"  "th_lyr"     {
/*capture:*/ svmat _`coef' 		// store average markdown for this bootstrap iteration in matrix
}
 
keep in 1/`B'

list _bl* _bm* _bk*


* Reshape and summarize kappa by year. Extract SD per year (or 5%, 95% CIs per year)? Plot this.
   
save ./data/temp/bootstraps, replace

sum m_kappa*, d
sum _kappa*, d

use ./data/temp/bootstraps, clear
sum _md_bb1, d 
list m_md_rts1 

sum _th_union*, d
 
 

** Standard errors and 5%-95% C.I.s

foreach var of varlist *_md*  *_mu* _th*   _b*  _scale*  *_kappa* *rho* {
egen se`var' = sd(`var')
egen p5`var' = pctile(`var'),p(5)
egen p95`var' = pctile(`var'),p(95)
}
 

global se_kappa_dq4 = se_kappa_dq4
global sem_kappa_dq4 = sem_kappa_dq4  

global se_kappa_rts_dq4 = se_kappa_rts_dq4
global sem_kappa_rts_dq4 = sem_kappa_rts_dq4  

 
di $se_kappa_dq4
di $sem_kappa_dq4

  
 *** Save SEs for tables

** Production function

 
** Store results to table

* markdown correlations 
foreach n in "a" "b" "c" {
forvalues m = 1/3 {
gen mdcorr = se_mdcorr`m'`n'
estpost su mdcorr 
est store mdcorr_`m'`n'_se 
drop mdcorr 
}
} 


* cobb-douglas
 
 local m = 1
foreach mod in "ols" "bb" "rts" {
rename (se_bl_`mod'`m' se_bk_`mod'`m'  se_bm_`mod'`m'  se_rho_`mod'`m'  se_thm_`mod'`m' 	 se_thl_`mod'`m'	se_scale_`mod'`m'    ) (bl bk bm rho thm   thl scale )

estpost su bl bm bk  rho
est store pf_`mod'`m'_se
estpost su thl thm scale  
est store th_`mod'`m'_se
drop bl bm bk    thm thl scale rho
}
 
 
** Markups & markdowns
 
sum se*
 
 forvalues m = 1/1 {
 	foreach mod in "bb" "ols" "rts" {
 foreach var in "mu" "md"   {
rename (se_`var'_`mod'`m' sem_`var'_`mod'`m' sea_`var'_`mod'`m' se_`var'p_`mod'`m' ) (av`var'  med`var'   ag`var' wedge`var')
}
estpost su medmd avmd medmu avmu 
est store md`m'_`mod'_se
drop medmd agmd avmd wedgemd medmu agmu avmu wedgemu 
}
}

 
 forvalues m = 1/1 {
foreach var in "p5" "p95" {
rename (`var'_md_bb`m'1  `var'_mu_bb`m'1 `var'm_md_bb`m'1  `var'm_mu_bb`m'1 `var'a_md_bb`m'1  `var'a_mu_bb`m'1  ) (avmd_bb`m'  avmu_bb`m' medmd_bb`m'  medmu_bb`m' agmd_bb`m'  agmu_bb`m'  )
estpost su avmd_bb`m' agmd_bb`m' medmd_bb`m' 
est store md`m'_bb_`var'
estpost su avmu_bb`m'  agmu_bb`m'  medmu_bb`m' 
est store mu`m'_bb_`var'
drop avmd_bb`m'  avmu_bb`m' medmd_bb`m'  medmu_bb`m' agmd_bb`m'  agmu_bb`m'
} 
}
  
* Time trends

rename (se_th_yr1 se_th_yr2 se_th_yr3 se_th_yr4 se_th_yr5 se_th_yr6 se_th_lyr ) ( th_yr1  th_yr2  th_yr3  th_yr4  th_yr5  th_yr6 th_lyr)

estpost su th_lyr
est store time1_se

estpost su   th_yr1  th_yr2 th_yr3  th_yr4 th_yr5 th_yr6     
est store time2_se


* Concentration

forvalues m = 1/2 {
rename (se_th_rail`m' se_th_tram`m' se_th_urban`m' se_th_nm1`m' se_th_nm2`m' se_th_nm3`m' ) ( th_rail  th_tram  th_urban  th_nm1 th_nm2 th_nm3) 
estpost su  th_rail   th_tram  th_urban  th_nm1 th_nm2 th_nm3 
est store conc`m'_se
drop th_rail  th_tram  th_urban  th_nm1 th_nm2 th_nm3 
}

* Collusion

sum se_th_car1 se_th_union1
 
forvalues m = 1/1 {
rename (se_th_car`m' se_th_union`m' ) ( th_car th_union) 
estpost su  th_car th_union 
est store col`m'_se
drop th_car th_union
}

forvalues m = 3/4 {
rename (  se_th_union`m' ) (  th_union) 
estpost su    th_union 
est store col`m'_se
drop   th_union
}

  
